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In this paper, we propose to use the ensemble Kalman smoother (EnKS) 
as linear least squares solver in the Gauss-Newton method for the large 
nonlinear least squares in incremental 4DVAR. The ensemble approach is 
naturally parallel over the ensemble members and no tangent or adjoint 
operators are needed. Further, adding a regularization term results in replacing 
the Gauss-Newton method, which may diverge, by the Levenberg-Marquardt 
method, which is known to be convergent. The regularization is implemented 
efficiently as an additional observation in the EnKS. Copyright (c) 2013 Royal 
Meteorological Society 
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1. Introduction 

4DVAR is a dominant data assimilation method used in 
weather forecasting centers worldwide. 4DVAR attempts 
to reconcile model and data variationally, by solving a 
very large weighted nonlinear least squares problem. The 
unknown is a vector of system states over discrete points 
in time, when the data are given. The objective function 
minimized is the sum of the squares of the differences 
of the initial state from a known background state at the 
initial time and the differences of the values of observation 
operator and the data at every given t ime point. In the 
weak-constraint 4DVAR (lTremolelll2007h . considered here, 
the model eiTor is accounted for by allowing the ending 
and starting state of the model at every given time point 
to be different, and adding to the objective function also 
the sums of the squares of those differences. The sums of 
the squares are weighted by the inverses of the appropriate 
error covariance matrices, and much of the work in the 
applications of 4DVAR goes into modeling those covariance 

matrices. 

In the incremental approach (ICourtier et al.\ [l994h . the 
nonlinear least squares problem is solved iteratively by 
using a succession of linear least square solutions. The 
major cost in 4DVAR iterations is in evaluating the 
model, tangent and adjoint operators, and solving large 
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linear least squares. A significant software development 
effort is needed for the additional code to implement 
the tangent and adjoint operators to the model and the 
observation operators. Strai ghtforward l inearization, called 
the incremental approach (ICourtier et a l. 1994), lea ds to 
the G a uss-Newton method for nonlinear least squares (iBelll 
119941; iTshimanga et al.\ |2008|) . However, Gauss-Newton 
iterations may not converge, not even locally. Finally, while 
the evaluation of the model operator is typically parallelized 
on modern computer architectures, there is a need to further 
parallelize the 4DVAR process itself. 

The Kalman filter is a sequential Bayesian estimation of 
the gaussian state of a linear system at a sequence of discrete 
time points. At each of the time points, the use of the 
Bayes theorem results in an update of the state, represented 
by its mean and covariance. The Kalman smoother simply 
considers all states at all time points from the beginning 
to be a large composite state. Consequently, the Kalman 
smoother is obtained from the Kalman filter by simply 
applying the same update as in the filter to the past states as 
well. Howe ver, historical l y, the f ocus was on efficient shor t 
recursions (iRauch efg/.l 119651: IStrang and Borrel |1997|) . 
similar the sequential Kalman filter. 

It is well known that weak constraint 4DVAR is 
equivalent to the Kalman smoother in the linear case. To 
apply the Kalman smoother in the nonlinear case, the 
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problem needs to be linearized, leading to variants of the 
extended Kalman filter and the Gauss-Newton method. Use 
of the Kalman smoother to solve the linear least squares in 
the Gauss-Newton method is known as the iterated Kalman 
smoother, and considerable improvements can be obtained 
against running th e Kalman smoother only once dBellll 19941 : 
iFisher et a/.ll2005l) . 

The Kalman filter and smoother require maintaining 
the covariance of the state, which is not feasible for 
large systems, such as in numerical weather prediction. 
Hence, the ensemble Kalm an filter (EnKF ) and ensemble 
Kalman smoother (EnKS) ( lEvensenl l2009l) use a Monte- 
Carlo approach for large systems, representing the state 
by an ensemble of simulations, and estimating the state 
covarian ce from the ensemb le. The implementation of the 
EnKS in lStroud et a/.l(l2010l) uses the adjoint model with the 
sh ort recursion s as in t h e KS. H owever, the implementations 
in iKhare et ai\ (|2008|) : lEvensen (.2009.) do not depend on 
the adjoint model and simply apply EnKF algorithms to the 
composite state over multiple time points. We use the latter 
approach here. 

The EnKF has become a competitive method for data 
assimilation. Consequently, combinations of ensemble and 
variational approaches have become of considerable recent 
interest. Estimating the background covariance for 4DVAR 
from an ensemble was one of the first connections 
dHamill and Snvden l2000t). and i t is now standard and 
became operational dWand l2010h . Gradient methods in 
the span of th e ensemble for one analysis cycle (i.e., 
3DVAR) include lZupanskil (l2005h . [Sakov et al\ (l20 1 2h (with 
squ are root EnKF as a l inear solver in Newton method), 
and iBocquet and Sakovl d2012h . who added regularization 
and use LETKF-like approach to minimize the nonlinear 
cost func t ion over line ar combinations of the ensemble. 
iLiu et al\ d2008l l2009h combine ensembles with (strong 
constraint) 4DVAR and minimize in the observation space. 
Their method, called Ens 4DVAR, doe s not n eed tangent 
or adjoint operators also. IZhang et al\ (l2009h use a two- 
way connection between EnKF and 4DVAR, to obtain the 
covariance for 4DVAR, and 4DVAR to feed the mean 
analysis into EnKF. However, ensemble methods for the 
solution of the 4DVAR nonlinear least squares problem 
itself or for the weak constraint 4DVAR do not seem to have 
been developed before. 

In this paper, we propose to use the ensemble Kalman 
smoother (EnKS) as linear least squares solver in 4DVAR. 
The ensemble approach is naturally parallel over the 
ensemble members. The rest of the computational work is 
relatively cheap compared to the ensemble of simulations, 
and parallel dense linear algebra libraries can be used. 
The proposed approach uses finite differences from 
the ensemble, and no tangent or adjoint operators are 
needed. To stabilize the method and assure convergence, 
a Tikhonov regularization term is added to the linear 
least squares, and the Gauss-Newton method becomes the 
Levelberg-Marquardt method. The Tiknonov regularization 
is implemented within EnKS as a computationa lly cheap 
additional observation (I Johns and Mandell l2008h . We call 
the resulting method EnKS-4DVAR. 

The paper is organized as follows. In Section|2] we review 
the formulation of 4DVAR. The EnKS for the incremental 
linearized squares problem is reviewed in Section |3] The 
new method without tangent operators is introduced in 
Section H] The modifications for the regularization and 
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the Levenberg-Marquardt method are presented in Section 
|5] Section |6] contains the results of our computational 
experiments, and SectionQis the conclusion. 

2. Incremental 4DVAR and the Gauss-Newton method 

We want to determine xq, . . . ,Xk, where Xi is the state at 
time i, from the background state, 



the model. 



Xo ~ Xh, 



Xi « Mi (x^-i) , 



and the observations 



y-i (xi) « Hi 



where Mi is the model operator, and Hi is the observation 
operator Quantifying the uncertainty by covariances, with 
Xo ~ Xh taken as {xq — x^) B~^ {xo — Xh) ~ 0, etc., we 
get the nonlinear least squares problem 

k 

i=l 
k 

EWVi - Hi (a;j)||^-i -^ min, (1) 



called weak-constraint 4DVAR dTremoletl2007l) . Originally 

in 4DVAR, Xi = Mi {xi^i); the weak constraint Xi w 
Mi (xi-i) accounts for model error 

The least squares problem ([T]i is solved iteratively by 
linearization. 

Ml {x,^i + Sx.^i) w M^ (x^-i) + M'i {xi-i)5xi-i, 
Hi {xi + 5xi) ^ Hi {xi) + H'i {xi) 6xi. 

For k vectors Ui,i~l...k, denote the composite vector 



UO:k 



Uo 



Uk 



In each iteration XQ-k ^- xo-k + Sxo;k, one solves the 
auxiliary linear least squares problem for the increments 

5xo.,k, 

\\xo + Sxo -aJblle-i 

fc 
+ X \\xi +5xi -7Wj(.Tj_i) - M'i (xi-i) Sxi-i\\q-i 



i=l 
k 



+ y^ hi - '^i {xi) - H'i [xi) 5xi\\p^-i -^ mill 

^ — ' » oxn-k 



(2) 



This is the Gauss -Newton method dBelll Il994t 
iTshimanga et al\ |2008|) for nonlinear sq uares, known 
in 4p VAR as the incremental approach (ICourtier et al.\ 
11994 . 
Denote 

ZO:k = Sxo.,k, Zb^Xb-Xo, 

mi= Mi{xi-i) - Xi, di^yi~Hi{xi) , (3) 
M, = M', (x,_i) , Hi = H', {xi) 
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and write the auxiliary linear least squares problem (|2) as and 



N aNT 






I I I N -I 



(10) 



l^o-^blls 

A- 



MjZi_i - mi 



iQ,- 



1 

2 



Hence, the analysis ensemble Z^. consists of linear 
combination of the forecast ensemble, which can be written 

N 



i=l 



mill 



fA\ as multiplying the ensemble by a transformation matrix T 



zN 



N rrN 



ryl\ r7i\ rri 



The function minimized in (|4]i is exactly the same as the 
one minimized in the Kalman smoother The Gauss-Newton 
method with the Kalman smoother as the linear least 
squares solver is known as the iterated Kalman smoother, 
and considerable improvements can be obtained against 
running the Kalma n smoother, applied to the l inearized 
problem, only once ( lBelllll994l:lFisher et a/.ll2005l) . 

3. Ensemble Kalman Filter and Smoother 

We prese nt the EnKF an d EnKS algorithms, essentially 
following lEvensenl (|2009|) . in a form needed to state our 
theorems. We start with a formulation of the EnKF in 
a notation suitable for extension to EnKS. The notation 
v^ ^ N (to, A) means that v^ is sampled from N (to, A) 
independently of anything else. The ensemble of states of 
the linearized model at time i, conditioned on data up to 
time j (that is, with the data up to time j already ingested), 
is denoted by 



where 



T^ = 1- 



N - 1 



Tf e M^><^, 



I - i^^ Af T 



1 



N aNT 



N ~ 1 



-A" A 



N 



Ri 



HiZ.^K^ -d, + w. 



i=l,N 



(11) 



(12) 



(13) 



with wf ^ N (di, Ri), and 

4^ = /f.Z,^_i(/-^) = [a,h...,af], 



-f^»^iK- 



1 ^ 



(14) 



j=i 






'«|j' ■ 



7^ 



'i\j 



where the ensemble member index £ always runs over i 
1, . . . ,N, and similarly for other ensembles. 



Remark 2 The matrix formula in the analysis step ^ is 
not efficient when the dimension of the data space is large. 
Using ( liOD and the Sherman-Morrisson-Woodbury formula 
({Hager 198^, we transform the inverse in into 



Algorithm 1 (EnKF) 1. Initialize 

2. For i = 1, . . . , fc, advance in time 
z-\i_i = M,zf_i|^_i + TO, + vl v\ 
followed by the analysis step 

■{Hizl^^_^-d^-w'i), wl^ 



[H,P^Hj + R.,y' =R-' 



..,N. 



(5) 



-A^ ( I -\ 



N -1 



N - 1 



(15) 



^iV(0,Q,), (6) 

R^)-' 

^7V(0,H,), (7) 



where 



jN 



1 



Using t[15\l requires only the solution of systems with the 
data error covariance matrix Ri (which is typically easy, 
and often Ri is diagonal) and of a system of the size N, the 
number of ensemble members. See \Mandel et al.l n200^} for 
details and operation counts. 

The EnKS is obtained by applying the same analysis step 
O as in the EnKF to the composite state Zo:i|i_ifrom time 
to i, conditioned on data up to time i — 1, 



-(^. 



N 



i TV — 1 '''*^"^ 

is the sample covariance, 

z^ 






ll^)(^. 



TV 
i\i — l 



-^1 a^)^ 



-ili-l-' 



(8) 



z, 



N 
0:?:|i-l 



yN 

z— 1 






yN ±_ 

^»|i-l^ 



in the place of Ziu^i. The observation term HiZ^, , — Di 



is the sample mean, and 1 is the vector of all ones size 
Nxl. 



becomes 



'i\i-l 



Remark 1 From (HJ), 
1 



,7V 



zN 



'K-ijV - 1 



11^ 



11^ 



Z. 



NT 

■i|i— 1' 

(9) 



[0, . . . , H,] .^0:,|i-l - A - -ffi-^iji^i - D 

Algorithm 2 (EnKS) Given Zb, 
1. Initialize 



■^o|o 



iV(zb,B) 



1,- 



,N. 



(16) 



(17) 
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2. For i — 1, . . . ,k, advance in time, 
followed by the analysis step 



(18) 






Po:i,0:lHo-AHo:lPo:i,0:lHQ.^ 



(19) 



(Jfo:^<,_i - A), A - N id,,Ri) , 



Algorithm 3 (Nonlinear EnKS) 1. Initialize 

a:g|o-7V(a;b,iB). 
2. For i = 1, . . . , fc, advance in time 

4\^-l-^M^[xl,^^^,)+v,, v,^NiO,Q,) (22) 
followed by the analysis step 



rN 



-TV 



-,W 



where i?o:i = [0, • . . , Hi], and P^-i oi is the sample 



-V-7V "v-iV /T-iiV /T-iA^ ^ TrpNxN 

^0:i\i - ^O-.ili-l-'- i ' J- I ^ ^ I 



covariance matrix of Z^^,-_y 

The following theorem allows a straightforward imple- 
mentation of the EnKS from the EnKF - the same trans- 
formation matrix is applied to the composite state from 
times to i, not just the last time i. Also, one can use 
a transformation matrix from another version of EnKF, 
such as the square root filter, e.g., LETKF (iHunt et al.\ 
l2007l) : lFertip et al\ ( |2007[) assume such relation a-priori for 
a related method based on LETKF. 



where 



Tf" ^I- 



N -1 



N -1 ' * 



11^ 

N 



A. 



NT 



Theorem 1 The EnKS satisfies 



yN _ yN rpN 

^0:i\i ~ ^OiilJ-l-'- i ■ 



(20) 



where T.^ is the transformation matrix Ml\ from the EnKF. 
Proof. We have 



,7V 



with the blocks 



TyN 



T}N 

^i,0 



pN 
^0,i 



jN 



,N 



N 



^(^- 



i-il )(^£k-i 



,T\T 



The terms in (|7]i become in ( fT9] l 



jN 



P" H P^ 



H 



pN „T 



n, 



y'-ili-lj 



A 



N 



[al,...,af] , 



Ri 

N 



1=1, N 



wf^N{0,R,), 
(23) 



a,- —Tiiix: 



H\i- 



i) - ^E^^ H.-i) ' 



i=i 



(24) 



andxi\i_i = X^^.^-^l/N. 

Using Theorem [T] it is easy to see that Algorithm |3] 
coincides with Algorithm|2]in the linear case, i.e., when Adi 
and "Hi are affine operators. 

4. Nonlinear EnKS-4DVAR method 

So far, the algorithm was relying on the linearized (i.e., 
tangent) model operators Mi and Hi and their adjoints. 

The linearized model Mi — M[ (xi-i) occurs only in 
advancing the time as action on the ensemble members 

5x^ — z^, 

MiZ-_^ + m, = M'i (xi-i) z-„i + Mi (xi^i) - Xi 

Approximating by finite differences based at Xi^i with step 
T > 0, we get 



^ Po:t.Q:iH = HiP^^H^ , 
and, similarly as in (|9|l, 

llT\ /, llT 



rN 



1 



pN __ 



N 



N 



Z, 



NT 



M^zf_^ 



Ml (xi_i + Tzf_i) - Ml [xi^i] 



(25) 



+ Ml [xi-i) - Xi 



(21) 

The result now follows by the comparison of (fT&t- (l2ll 
with dTt-lfTTT). □ 

When the original, nonlinear operators instead of the 
linearizations are used, we get the nonlinear EnKS method, 
which is common and useful in practice, even if it may not 
be justified theoretically. This method is obtained from the 
linear EnKS by replacing ( fTSl l and (fT4l i by their original, 
nonlinear versions. It operates on the original ensemble 

of the states X^ = [x^] . rather than on the increments 

z^ = 5x^. 



Thus, advancing the linarized model in time requires iV + 1 
evaluations of Mi, at Xi^i and Xi^i + t5x'^_i. 

The observation matrix Hi occurs only in the action on 
the ensemble, 

HiZ«=[Hizl...,HiZ^]. 

Approximating by finite differences based at Xi, with step 
T > 0, we have 



Hz 



Hi (Xi-l + Tzf) ~ Hi (Xi^l) 



(26) 



Copyright © 2013 Royal Meteorological Society 
Prepared using qjrms4,cls 



Q. J. R. Meteorol. Soc. 00:[T](9](2013) 



4DVAR by ensemble Kalman smoother 



Thus, evaluating the action of the linarized observation 

requires iV + 1 evaluations of Hi, at Xi_i and Xi-i + 
t 

T-4-1- 

Here is the overall method. First, initialize 

xq = Xb, Xi = Mi {xi-i) , i = 1, . . . , fc, 

if not given already. One iteration (|2) of the incremental 
4DVAR is then implemented as follows. 



■ ,Xk: 



Algorithm 4 (EnKS-4DVAR) Given xq, 

1. Initialize Zqiq ~ A^ (zb, B) following 0, with Zh = 0. 

2. For i = 1, . . . ,k, advance z^in time following M8\ 
with the linearized operator approximated from ( 1251 ), 



M.. 



i\i-l 



ix^-i + Tzl_^y^_A - Mt {xi^i) 



(27) 



+ M^{x,.i)-x, + vl wf - TV (0, Q,) , 

followed by the analysis step i20\l with the transformation 

matrix T^ computed from ( 1721 ) and with the matrix-vector 

products HiZi approximated from { 

3. Update 



Xi 



1 ^ 
N ^-^ 



N 



■j|fc' 



0, 



Note that for small r, the resulting method is asymptoti- 
cally equivalent to the method with the derivatives. Amaz- 
ingly, it turns out that in the case when t = 1, we recover the 
standard EnKS applied directly to the nonlinear problems, 
that is, with the linearized advance in time (|6]l replaced 
by application of the original, nonlinear operator A^^. In 
particular, the incremental 4DVAR does not converge unless 
it is already at a stationary point, because each iteration 
delivers the same result, up to the randomness of the EnKS. 

Theorem 2 If r = 1, then one step of EnKS-4DVAR 
(Algorithm^ is exactly the nonlinear EnKS (Algorithm\^. 
In particular, the result of the step does not depend on the 
previous iterate. 



Proof. Indeed, ( |27] | becomes 



Af. 



+ z.. 



M^ 



-i]!-! 



1 



hence 



Xi + z. 



-A(» 



/x. 



■V, 



i-l) X._ 
Xi-1 + ^i_i|i_i 



M^ 



(28) 



which is exactly the same as advancing the ensemble 
member £ following ( |22] ) with : 



Similarly, (fl4l i becomes with t = 1 

H^ (x, + zf|^_ J - n, (Xi) 
"'= 1 

1 ^ m{x. + zi\.,_^-H.ix^ 
N ^ 1 



+ z: 



i-l|i-l- 



J = l 



Hi Xj + zf 



1 ^ 



N 
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(30) 



which is exactly the same as ( l24l i with xi^_^ = Xi + zi^_-^. 
Finally, ( fT3] ) becomes using (O, 



-ffiZi|i_i - di 

Ui (xi + ^f|j_i j - Tit {xi) 
^ 1 



-[y^-n,ix,)] (31) 
(32) 



which is exactly the same as in (1231 ). □ 

5. Tikhonov regularization and the 
Levenberg-Marquardt method 

The Gauss-Newton method may diverge, but convergence 
to a stationary point of ([T]i can be recovered by a 
control of the step Sx. Adding a constraint of the form 
WSxiW <£ leads to gl obally convergent trust region methods 
(iGratton et al.\ l2013l) . Here, we add Sxi in a Tikhonov 
regularization term of the form 7 ||(5xi||g-i, which controls 
the step size as well as rotates the step direction towards 
the steepest descent, and obtain the Levenberg-Marquardt 
method xo;k <— xo-.k + ^xo-.k, where 



\\Sxo - Zblls-i + X! 1'^^* ~ Mi6xi-i - m,||Q-i 

4 = 1 



+ y^ \\dt - HiSxiWjf^ 

4=1 

fe 
+ j^\\Sx, 



i=0 



1 — !> mm 



(33) 



Under suitable technical assumptions, the Levenberg- 
Marquardt method is guaranteed to converge globally 
if the regularization parameter 7 > is large enough 
dOill and Murravlll978t rOsbomelll976l) . Estimates for the 
convergence of the Levenberg-Marquardt method in the 
case when the linear system is solved only approximately 
exist dWright and Holj[l985h. 

Similarly as in Ijohns and Mandell (l2008h . we interpret 
the regularization terms 7 |j(5a;i|jg-i in ( l33T l as arising from 

additional independent observations Sxi ^ A(0,7~^S'i) 
Because the additional regularization observations Sxi « 
are independent of the other observations and the state, 
separately, resulting in the mathematically equivalent but 
often more efficient two-stage method - simply run the 
EnKF analysis (|7]) twice, and apply both transformation 
matrices in turn following (l20t . With the choice of Si 
as identity or, more generally a diagona l matrix, the 
implem entation following ( fTSl l is efficient; see lMandel et al.\ 
(|2009|) for operation counts. 

Note that unlike in iJohns and Mandell (|2008|) . where the 
regularization was applied to a nonlinear problem and thus 
the sequential data assimilation was only approximate, here 
the EnKS is run on the auxiliar linearized problem (l3Jt , so 
all distributions are gaussian and the equivalence of solving 
( l33b at once and assimilating the observations sequentially 
is statistically exact. 
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Figure 1. The Lorenz attractor, initial values x(0) = 1, j/(0) = 1, 
z(0) = 1, time step dt = 0.1. 



6. Computational results 

6.1. Lorenz 63 model 

We first show an example without model error, where 

convergence is achieved with 7 = 0(algo|4l). 

We consider the Lorenz 64 equations (iLorenzl Il963l) . 
a simple dynamical model with chaotic behaviour The 
Lorenz equations are given by the nonlinear system 



dx 

'dt 
dy 

dz 
It 



a{x - y) 



xy — f3z 



where x = x{t), y = y{t), z = z{t) and a, p, (3 are 
parameters, which in these experiments are chosen to have 
the values 10, 28 and 8/3 respectively. The system is 
discretized using the fourth-order Runge-Kutta method. In 
([B, we choose 




Figure 2. Nonlinearity of the Lorenz 63 model. The values of x (i + 
1), y{t + 1) and z{t + 1) change quickly as a function of x{t) = 1, 
y{t) = 1, and varying z{t). 



Iteration 


1 


2 


3 


4 


5 


6 


RMSE 


20.16 


15.37 


3.73 


2.53 


0.09 


0.09 



Table 1 . Norm of the root mean square error in Gauss-Newton iterations 
with EnKS as linear solver. 





B ^ al diai 



g 1, 



1 1 
4' 9 



R, = </, 



'Hi{x,y,z) 



ix ,y ,z ) 



In the experiments below, we assume perfect model, 
therefore Q^ ~ el, e w 0. We put 0-5 = 1, cr,. = 1, and e = 
0.0001. 

As can be seen in Figures [1] and [2] the Lorenz attractor is 
fully nonlinear It has two lobes connected near the origin, 
and the trajectories of the system in this saddle region 
are particularly sensitive to perturbations. Hence, slight 
perturbations can alter the subsequent path from one lobe 
to the other In Figure |2] we keep x{t) and y{t) constant 
and we vary just z{t), then we compute the state at time 
t -\- \, this figure shows the non linear dependence between 
the different components of the state at time t -f 1 and the 
third component of the state at time t. 

To evaluate the performance of the method, we use the 
twin experiment technique. That is, an integration of the 
model is chosen as the true state. We then obtain the 
data yi by applying the observation operator T-Li to the 
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Figure 3. The three components x, y, z of the truth and five iterations of 
EnKS-4DVAR. The initial conditions for the truth are a::(0) = 1,3/(0) = 1, 
and z{0) = 1, time step dt = 0.1, observations are the full state at each 
time, ensemble size is 100. 
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Figure 4. Root mean square error of EnKS-4DVAR iterations. Tlie 
problem setting is the same as in Fig.|3] 



Figure 5. Evolution of the 40-variable Lorenz-96 system at site 20. The 
red hne is unperturbed forecast. The black lines are an ensemble of 50 
forecasts, which start from slightly perturbed initial conditions. 



truth and then adding a gaussian perturbation A^(0, Ri). 
Similarly, the background Xb is sampled from the gaussian 
distribution with the mean equal to the initial conditions and 
the covariance B. Then we try to recover the truth using the 
observations yi and the background Xf,- 

Figure |3] reports simulation results for assimilating 
observations over 50 assimilation cycles, using the Hybrid 
4DVAR and nonlinear EnKS method. Cycles are separated 
by a time interval of dt = 0.1. Figure |4] and shows the root 
mean square error (RMSE) between the sample posterior 
mean and the true state of the system. 

As can be seen from Table [T] five iterations were enough 
for the method to converge. Note that the error does 
not converge to zero, because of the approximation and 
variability inherent in the ensemble approach. 

6.2. Lorenz 96 model 

We now consider the effect of the model error in the 

algorithmic 

The Lorenz 96 model dLorenzl [20061) is defined by the 
system of differential equations 



dxj_ 
dt 



1 



K 



{xj~i{x. 



J + 1 



Cj-2 



F), 



J = 

Xq 



1, . . . , 40, with cyclic boundary conditions a;_i = 0:39, 
= XiQ, xai = xi. This model behaves chaotically in 



the case of external forcing F = 8. The first term of 
right-hand side simulates advection, and this model can 
be regarded as the time evolution of a one-dimensional 
quantity on a constant latitude circle, that is, the subscript 
corresponds to longitude. The PDF is discretized using a 
fourth-order Runge-Kutta method. For the results below 
K = 1, F = 8 and dt — 0.01. Figure |5] shows the chaotic 
dynamics of the Lorenz 96 system. For the tests we 
took the parameters B = (T^diag{l, ..., -^, ...), Ri = tr^/. 



(TqC, where Cij 



Xf, 7W. 



is the Lorenz 96 model, Qi = 



exp{-\z-jm,L 



dt 



, o-fc 



l,cr.r 



1, (7q = 0.1, and the ensemble size N = 40. We use again 
the twin experiments technique, the true state is equal to 
an integration of the model plus a gaussian perturbation 




5 10 t5 20 25 30 35 40 45 50 



Figure 6. EnKS-4DVAR for the Lorenz 96 system. The truth at t = is 
generated randomly. The figure shows the truth and the filter by EnKS- 
4DVAR at xio{t), X2oit), and X3o(t). 



A(0, Qt). Figures |6] and Q illustrate EnKS-4DVAR on this 
problem. 

6.3. Example where Gauss-Newton does not converge 

We now show that the algorithm |4] may not be convergent 
and that Tikhonov regularization may be needed in some 
circumstances. The following academic example illustrates 
this fact. 

The Gauss Newton method for nonlinear least squares 
is not globally convergent, but convergence to a stationary 
point of any least square problem can be recovered 
by using the Levenberg-Marquart control. Consider the 
following example, which requires Levenberg-Marquart 
regularization to converge. The objective function to 
minimize is 

J{xo,xi) = {XQ - 2)2 + (3 + x\f + -(xo - xif (34) 
where q is small, q = 0.000001. 
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Figure 7. Root mean square error between the truth and EnKS-4DVAR for 
the Lorenz 96 system and the same setup as in Figure |6] 
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iteration number 
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Figure 8. The RMSE between Gauss-Newton iterations and local 
minimum {xq,x'^) = (0.419, 0.419)(7 = 0, top), and the root mean 
squared eiTor between the Levenbert-Marquardt iteration and the local 
minimum (xq, Xj)(7 = 200, bottom). 



This problem could be seen as a 4DVAR problem where 
the state at time is xq, the background state is Xb = 2, 
the background covariance B = I, there is one time step, 
the state at time t ~ 1 is xi, the model A/i = /, and the 
model is perfect. Qi = 0.000001 « 0, observation operator 
7ii{x) = ~x^ and observation error covariance is i?i = /. 

Figure [8] shows the iterations of the EnKS-4DVAR 
method applied to the problem |34] seen as 4Dvar problem, 
in two cases: when 7 = ( Gauss Newton iteration) the 
method does not converge, and for 7 = 200 ( Levenberg- 
Marquart iteration), the method seems to converge to the 
local minimum (a;^, x\) ^ (0.419, 0.419). 

7. Conclusion 

The EnKS-4DVAR method was formulated and shown to 
be capable of handling strongly nonlinear problems, and it 
converges to the nonlinear least squares solution in a small 
number of iterations. Its performance on large and realistic 
problems will be studied elsewhere. 
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